tHIS IS THE SIMULATION i HAD IN THE stc PROJECT. BUT TI GIVES poor identification of parameters, with large sets. Even when using const',HC3orNeweyWest` estimators.
Simulation data comes from C:\Users\Emi\OneDrive\UWA PhD\MonPolicy StockMkt\Results\Simulations\NoCommShockZ\CalCrisis-HARRVBig-NoZ-1016-Tres\simVARs.RDS[[1]]
see details in that folder.
INIT
Code
#required packages to make a packagepacman::p_load(devtools, roxygen2, testthat)#dependecies#pacman::p_load(tidyverse, data.table, DT, plotly, gridExtra, latex2exp, here, logr, future.apply, furrr)#pacman::p_load(dplyr,parallel,sandwich,here,future,logr,profvis,lobstr,parallel)#https://ourcodingclub.github.io/tutorials/writing-r-package/#create_package#https://combine-australia.github.io/r-pkg-dev/functions.html#load_all()#run these two lines and then library(GridMaker) in normal rstudio to install package# load_all loads a package. It roughly simulates what happens when a package is installed and loaded with library().devtools::load_all()devtools::install()
The simulated model is a bivariate SVAR(2) as follows \[
\begin{equation*}
\begin{bmatrix}
1 & -\beta \\
-\alpha & 1
\end{bmatrix}
\begin{bmatrix}
i_t \\
s_t
\end{bmatrix}
=
\begin{bmatrix}
0.3 & 0.3 \\
-0.2 & 0.6
\end{bmatrix}
\begin{bmatrix}
i_{t-1} \\
s_{t-1}
\end{bmatrix}
+
\begin{bmatrix}
-0.2 & 0.1 \\
0.1 & 0.1
\end{bmatrix}
\begin{bmatrix}
i_{t-2} \\
s_{t-2}
\end{bmatrix}
+
\begin{bmatrix}
\varepsilon_t \\
\eta_t
\end{bmatrix}
\end{equation*}
\] where the impact coefficients are set similar to our findings for the before-ELB sample: \(\alpha=-5,\;\beta \times 100 = 1\) . The error terms \(\left\{ \varepsilon_{t}, \eta_{t}, \right\} _{t=1}^{T}\) are independently normally distributed with mean zero and time-varying variances \(\sigma^{2}_{\varepsilon; t}\), \(\sigma^{2}_{\eta; t}\).
1. ESTIMATION VAR
Code
varModel <-VAR(cbind('i'= sam$i, 's'= sam$s), p =2, type ="const")Res =residuals(varModel)
Code
# Use for computing the moment functionOmega <-cbind('ii'= Res[, 'i']*Res[, 'i'], 'is'= Res[, 'i']*Res[, 's'],'ss'= Res[, 's']*Res[, 's'])Y <-cbind( 'omega_ii'= Omega[, 'ii'], 'omega_ss'= Omega[,'ss'], 'omega_is'= Omega[,'is'] )# Compute the moment function at each t. # INPUT: a named vector with values for the structural parameters# OUTPUT: a TxG matrix where each row corresponds to the value of the moment function on that obseration (date)# an each column the value fot ehfat moment equationGET_MOM_VALUES <-function(t0, ...) {list2env(as.list(t0),envir =environment())#Y <- cbind( 'omega_ii' = Res[,'i']^2, 'omega_ss' = Res[,'s']^2, 'omega_is' = Res[,'i'] * Res[,'s'] )#Y %*% c(-theta0['alpha'],# -theta0['beta'],# 1 + theta0['alpha'] * theta0['beta'])#fT as.matrix(-alpha*Omega[, 'ii'] - beta*Omega[, 'ss'] + (1+alpha*beta) * (Omega[, 'is']))}
#### Search settings ----NCORES = params$NCORESSTEPSPERCORE =200# Size of adjacent cells to explores is NCORES * STEPSPERCORE# PARAMS_CFG a tibble where each row defines properties of the parameters in the name column# stepSize: size of the step for each parameter. You can have different sizes (and precision) for each paramter.# minVal and maxVal: Limits of the parameter space for each parameters. A list where each element defines the minimum and maximum for each parameter# ORDER: the order of the parameter has to match the order in theta (as used by the get_b() fcn)PARAMS_CFG <-tribble(~name, ~iniVal, ~stepSize, ~minVal, ~maxVal,'alpha', 0.3, 1, -20, 20,'beta', -0.2, 0.005, -0.03, 0.03) %>%mutate('N_TICKS'= ((maxVal-minVal)/stepSize)+1)
Code
NEQS =1MODELEQS <-list(list('Y'= Y, 'get_b'=NULL))#'X'=X, 'Z'=Z))NCORES = params$NCORES# Level of condfidence of the S set. The closer to 1, the more robust the search is to local minima but it takes longer.SPVALUE =0.95# R is NEQS*ncol(Z) where ncol(Z) is the number of instruments# By default this matrix has to be automatically built depending on the dimensions of MODELEQS.# The user can alternatively provide this matrixR =1if (exists('Z', mode='numeric')) { R <-rbind(cbind(matrix(1, nrow=ncol(Z), ncol=ncol(Z))) )}theta0 = PARAMS_CFG$iniValnames(theta0) = PARAMS_CFG$namehead(GET_MOM_VALUES(theta0))
expResult <-Stage2_rover_Perseverance(iniVal=PARAMS_CFG$iniVal, typeVarF ="NeweyWest", WTLags=1,HACprewhite =TRUE, localSearch =FALSE) # updates the Grid by in place
[1] "Local search? FALSE"
[1] "12:44. Do a GLOBAL search. nextRows contains all row numbers in the Grid."
[1] "Calling garbage collector before repeat to start exploration"
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 1258688 67.3 2345992 125.3 2345992 125.3
Vcells 3297808 25.2 8388608 64.0 5685111 43.4
[1] "12:44(top of repeat loop). Starting an exploration"
[1] "12:44hs. (inside for loop). Exploring 518 points. The first point is [-19.7, -0.03]. The last point is [19.3, 0.025]."
[1] "12:44hs. (inside for loop). Finished computing test values for this batch. Now about to update the Grid"
[1] "12:44hs. (inside for loop). Saving explored rows to file."
[1] "Calling garbage collector after finishing this batch in the\n for loop "
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 1345668 71.9 2345992 125.3 2345992 125.3
Vcells 3497108 26.7 8388608 64.0 5685111 43.4
# A tibble: 2 × 4
time length_minutes npoints points_per_minute
<dttm> <dbl> <dbl> <dbl>
1 2026-01-22 12:44:31 NA NA NA
2 2026-01-22 12:44:47 0.263 518 1968.
[1] "Perseverance has landed."
Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
Warning in geom_point(aes(x = -5, y = 0.01), size = 3, colour = "black", : All aesthetics have length 1, but the data has 155 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
a single row.
Code
pGplot
Warning in geom_point(aes(x = -5, y = 0.01), size = 3, colour = "black", : All aesthetics have length 1, but the data has 155 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
a single row.
Source Code
---title: "Unidentified case"author: "Emi Carlevaro"date: todaydate-meta: 14/01/2026editor: mode: sourceformat: html: code-fold: true code-tools: true self-contained: true page-layout: full toc: trueparams: IDSAMPLE: "" NCORES: 6 LOCALSEARCH: FALSE---tHIS IS THE SIMULATION i HAD IN THE stc PROJECT. BUT TI GIVES poor identification of parameters, with large sets.Even when using `const', `HC3` or `NeweyWest` estimators.Simulation data comes from `C:\Users\Emi\OneDrive\UWA PhD\MonPolicy StockMkt\Results\Simulations\NoCommShockZ\CalCrisis-HARRVBig-NoZ-1016-Tres\simVARs.RDS[[1]]`see details in that folder.## INIT```{r, eval=FALSE}#required packages to make a packagepacman::p_load(devtools, roxygen2, testthat)#dependecies#pacman::p_load(tidyverse, data.table, DT, plotly, gridExtra, latex2exp, here, logr, future.apply, furrr)#pacman::p_load(dplyr,parallel,sandwich,here,future,logr,profvis,lobstr,parallel)#https://ourcodingclub.github.io/tutorials/writing-r-package/#create_package#https://combine-australia.github.io/r-pkg-dev/functions.html#load_all()#run these two lines and then library(GridMaker) in normal rstudio to install package# load_all loads a package. It roughly simulates what happens when a package is installed and loaded with library().devtools::load_all()devtools::install()``````{r}#| warning: false#| message: false#| error: false#| include: true#| echo: truelibrary(vars)library(tidyverse)library(plotly)library(iStabi)``````{r}sam <- iStabi::data_sim_biVAR```## 0. THE MODELThe simulated model is a bivariate SVAR(2) as follows$$\begin{equation*} \begin{bmatrix} 1 & -\beta \\ -\alpha & 1 \end{bmatrix} \begin{bmatrix} i_t \\ s_t \end{bmatrix} = \begin{bmatrix} 0.3 & 0.3 \\ -0.2 & 0.6 \end{bmatrix} \begin{bmatrix} i_{t-1} \\ s_{t-1} \end{bmatrix} + \begin{bmatrix} -0.2 & 0.1 \\ 0.1 & 0.1 \end{bmatrix} \begin{bmatrix} i_{t-2} \\ s_{t-2} \end{bmatrix} + \begin{bmatrix} \varepsilon_t \\ \eta_t \end{bmatrix}\end{equation*}$$where the impact coefficients are set similar to our findings for the before-ELB sample: $\alpha=-5,\;\beta \times 100 = 1$ . The error terms $\left\{ \varepsilon_{t}, \eta_{t}, \right\} _{t=1}^{T}$ are independently normally distributed with mean zero and time-varying variances $\sigma^{2}_{\varepsilon; t}$, $\sigma^{2}_{\eta; t}$.## 1. ESTIMATION VAR```{r}varModel <-VAR(cbind('i'= sam$i, 's'= sam$s), p =2, type ="const")Res =residuals(varModel)``````{r}# Use for computing the moment functionOmega <-cbind('ii'= Res[, 'i']*Res[, 'i'], 'is'= Res[, 'i']*Res[, 's'],'ss'= Res[, 's']*Res[, 's'])Y <-cbind( 'omega_ii'= Omega[, 'ii'], 'omega_ss'= Omega[,'ss'], 'omega_is'= Omega[,'is'] )# Compute the moment function at each t. # INPUT: a named vector with values for the structural parameters# OUTPUT: a TxG matrix where each row corresponds to the value of the moment function on that obseration (date)# an each column the value fot ehfat moment equationGET_MOM_VALUES <-function(t0, ...) {list2env(as.list(t0),envir =environment())#Y <- cbind( 'omega_ii' = Res[,'i']^2, 'omega_ss' = Res[,'s']^2, 'omega_is' = Res[,'i'] * Res[,'s'] )#Y %*% c(-theta0['alpha'],# -theta0['beta'],# 1 + theta0['alpha'] * theta0['beta'])#fT as.matrix(-alpha*Omega[, 'ii'] - beta*Omega[, 'ss'] + (1+alpha*beta) * (Omega[, 'is']))} ``````{r, eval=FALSE}GET_B <- function(theta0) { c(-theta0['alpha'], -theta0['beta'], 1 + theta0['alpha'] * theta0['beta'])}```## 2. SEARCH### Search settings```{r}#### Search settings ----NCORES = params$NCORESSTEPSPERCORE =200# Size of adjacent cells to explores is NCORES * STEPSPERCORE# PARAMS_CFG a tibble where each row defines properties of the parameters in the name column# stepSize: size of the step for each parameter. You can have different sizes (and precision) for each paramter.# minVal and maxVal: Limits of the parameter space for each parameters. A list where each element defines the minimum and maximum for each parameter# ORDER: the order of the parameter has to match the order in theta (as used by the get_b() fcn)PARAMS_CFG <-tribble(~name, ~iniVal, ~stepSize, ~minVal, ~maxVal,'alpha', 0.3, 1, -20, 20,'beta', -0.2, 0.005, -0.03, 0.03) %>%mutate('N_TICKS'= ((maxVal-minVal)/stepSize)+1)``````{r}NEQS =1MODELEQS <-list(list('Y'= Y, 'get_b'=NULL))#'X'=X, 'Z'=Z))NCORES = params$NCORES# Level of condfidence of the S set. The closer to 1, the more robust the search is to local minima but it takes longer.SPVALUE =0.95# R is NEQS*ncol(Z) where ncol(Z) is the number of instruments# By default this matrix has to be automatically built depending on the dimensions of MODELEQS.# The user can alternatively provide this matrixR =1if (exists('Z', mode='numeric')) { R <-rbind(cbind(matrix(1, nrow=ncol(Z), ncol=ncol(Z))) )}theta0 = PARAMS_CFG$iniValnames(theta0) = PARAMS_CFG$namehead(GET_MOM_VALUES(theta0))``````{r}verify_global(PARAMS_CFG, Y)Times =NROW(Res) # time periodsgenerate_global(PARAMS_CFG, SPVALUE)```### Global search grid```{r}gridLines <-build_lines(PARAMS_CFG, 'iniVal')Grid <-build_grid(PARAMS_CFG, gridLines)```### Start search```{r}expResult <-Stage2_rover_Perseverance(iniVal=PARAMS_CFG$iniVal, typeVarF ="NeweyWest", WTLags=1,HACprewhite =TRUE, localSearch =FALSE) # updates the Grid by in place#saveRDS(Grid, 'Grid_Unidentified_model_WTLagsMinus1.rds')``````{r}head(Grid)``````{r}Grid[order(qLLStab)][1:10]```## 3. BUILD SETS### Critical valuesCritical values for $qLL-\tilde{S}$ for `r KZ` moment condition from table VII Suplementary material, p3 Magnusson Mavroeidis (2014).THe S-test follows a Chi-Square distribution with **`r KZMKX`** degrees of freedom (equatl to the number of moment conditions).`Grid` has `r NROW(Grid)` rows.The number of *fixed* identified parameters is **`r NFIXPARAMS`**.The number of *estimated* parameters under the null is **`r NSIPARAMS`**.```{r CVs}NSIPARAMS = 0 # not needed in the new version of the packagePCONFS = c(0.80, 0.85)get_CVs <- function() { CVs <- list('qLLStab'=NA, 'S'=NA, 'genS_qLL'=NA) CVs$qLLStab <- filter(iStabi::data_CVs$qLLStab, DG_F == KZ & P_CONF %in% PCONFS) CVs$genS_qLL <- filter(iStabi::data_CVs$genS_qLL, MOMENT_CONDITIONS == KZ & STRONGLY_IDENTIFIED_PARAMETERS == NSIPARAMS & P_CONF %in% PCONFS) CVs$S <- tibble('P_CONF' = PCONFS, 'CUTOFF_TOP' = map_dbl(PCONFS, ~qchisq(.x, KZMKX))) CVs}CVs <- get_CVs()```Critical values are:::: {.panel-tabset}### S```{r}CVs$S```### qLLStab```{r}CVs$qLLStab```### genS_qLL```{r}CVs$genS_qLL```:::### Data for plotting```{r}build_set <-function(setName) { thisCV <- CVs[[setName]] statMax =filter(thisCV, P_CONF == PCONFS[2])$CUTOFF_TOP thisPconf = thisCV$P_CONF thisCV = thisCV$CUTOFF_TOP theseCols <-c(PARAMS_CFG$name, setName) aSet = Grid[get(setName) <= statMax, ..theseCols][, sigAt := thisPconf[ 1+findInterval(get(setName), thisCV,left.open=TRUE)]] aSet}``````{r}Sets <-list('qLLStab'=build_set('qLLStab'),'S'=build_set('S'),'genS_qLL'=build_set('genS_qLL'))saveRDS(Sets, 'Sets_Unidentified_model.rds')```## 4. PLOTTING### Grid for plottingForce the starting values to be a multiple of the step size. This avoids problems and is clear what the true parameter space is.### PlotsTrue parameter vector $ \color{blue}{\theta} =(\alpha = -5, \beta = 0.01)$#### Global searchScatterplot for all searched values (useful for the global search)##### S plot```{r}# INTERACTIVE PLOTs ----#### S-plot ----pSbase <-plot_ly(data = Sets$S[sigAt <=0.80, ],x =~alpha,y =~beta,type ="scatter",mode ="markers",marker =list(color ="grey")) %>%layout(title ="S-set (90% confidence)",xaxis =list(title = plotly::TeX("\\alpha"),range =c(-20, 20),fixedrange =TRUE ),yaxis =list(title = plotly::TeX("\\beta"),tickvals =seq(from =-0.03, to =0.03, by =0.005),range =c(-0.02, 0.02),fixedrange =TRUE ) ) %>%add_trace(x =c(-5),y =c(0.01),type ="scatter",mode ="markers",marker =list(symbol ="square", size =12, color ="red"),inherit =FALSE,showlegend =FALSE ) %>%config(mathjax ="cdn")pSbase %>%add_trace(type='scatter', mode='markers', marker =list(color='blue'))```##### qLL-Stab plot```{r}# INTERACTIVE PLOTs ----#### S-plot ----pSbase <-plot_ly(data = Sets$qLLStab[sigAt<=0.80,], x=~alpha, y=~beta,color='grey')%>%layout(title ="qLL-Stab (90% confidence)",xaxis =list(title = plotly::TeX("\\alpha"),range =c(-20, 20),fixedrange =TRUE ),yaxis =list(title = plotly::TeX("\\beta"),tickvals =seq(from =-0.03, to =0.03, by =0.005),range =c(-0.02, 0.02),fixedrange =TRUE ) ) %>%add_trace(x =c(-5),y =c(0.01),type ="scatter",mode ="markers",marker =list(symbol ="square", size =12, color ="red"),inherit =FALSE,showlegend =FALSE ) %>%config(mathjax ="cdn")pSbase %>%add_trace(type='scatter', mode='markers', marker =list(color='blue'))```#### genS-qLL plot```{r}# INTERACTIVE PLOTs ----pSbase <-plot_ly(data = Sets$genS_qLL[sigAt<=0.80,], x=~alpha, y=~beta,color='grey')%>%layout(title ="genS_qLL (80% confidence)",xaxis =list(title = plotly::TeX("\\alpha"),range =c(-20, 20),fixedrange =TRUE ),yaxis =list(title = plotly::TeX("\\beta"),tickvals =seq(from =-0.03, to =0.03, by =0.005),range =c(-0.02, 0.02),fixedrange =TRUE ) ) %>%add_trace(x =c(-5),y =c(0.01),type ="scatter",mode ="markers",marker =list(symbol ="square", size =12, color ="red"),inherit =FALSE,showlegend =FALSE ) %>%config(mathjax ="cdn")pSbase %>%add_trace(type='scatter', mode='markers', marker =list(color='blue'))```#### Static plot: qLL-Stab```{r}pBase <-ggplot(data=Sets$qLLStab[sigAt<=0.80,], mapping=aes(x=alpha, y=beta)) +scale_x_continuous(name=TeX(r"(\\alpha)"), limits=c(-20, 20), breaks=seq(from=-20, to=20, by=2)) +scale_y_continuous(name=TeX(r'($\\beta$)'), limits=c(-0.03, 0.03), breaks=seq(from=-0.03, to=0.03, by=0.01)) +theme_bw() +theme(panel.grid=element_line(colour='#999999', linetype='14'), panel.grid.minor=element_line(colour='white'), text=element_text(size=10), axis.title.x =element_text(colour='red'),axis.title.y =element_text(colour='blue', size=rel(0.9), hjust=0.9),axis.text.x =element_text(size=rel(1.25), angle=90, vjust=0.5), axis.text.y =element_text(size=rel(1.25))) +geom_vline(xintercept=0, colour ='black', linetype ='solid') +geom_hline(yintercept=0, colour ='black', linetype ='solid') +# Highlightinh the TRUE valuegeom_point(aes(x=-5, y=0.01), size=3, colour='black', fill='yellow')# ScatterplotpGplot <- pBase +geom_point(aes(colour=qLLStab))ggsave('Unidentified_model_qLLStab_set.png', plot=pGplot, width=6, height=5)pGplot```